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Abstract 

Reaction-diffusion models are used to describe systems in fields as diverse as 
physics, chemistry, ecology and biology. The fundamental quantities in such 
models are individual entities such as atoms and molecules, bacteria, cells 
or animals, which move and/or react in a stochastic manner. If the number 
of entities is large, accounting for each individual is inefficient, and often 
partial differential equation (PDF) models are used in which the stochastic 
behaviour of individuals is replaced by a description of the averaged, or 
mean behaviour of the system. In some situations the number of individuals 
is large in certain regions and small in others. In such cases, a stochastic 
model may be inefficient in one region, and a PDF model inaccurate in 
another. To overcome this problem, we develop a scheme which couples a 
stochastic reaction-diffusion system in one part of the domain with its mean 
held analogue, i.e. a discretised PDF model, in the other part of the domain. 
The interface in between the two domains occupies exactly one lattice site 
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and is chosen such that the mean held description is still accurate there. This 
way errors due to the hux between the domains are small. Our scheme can 
account for multiple dynamic interfaces separating multiple stochastic and 
deterministic domains, and the coupling between the domains conserves the 
total number of particles. The method preserves stochastic features such as 
extinction not observable in the mean held description, and is signihcantly 
faster to simulate on a computer than the pure stochastic model. 

Keywords: Reaction-dihusion system. Stochastic model. Hybrid model, 
Fisher-Kolmogorov equation, Lotka-Volterra equation 
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1. Introduction 


Random effects due to finite numbers of components are ubiquitous in 
reaction-diffusion systems. Within this context, much research has been 
done, for instance, on the robustness and response to noise in gene regulatory 
networks. The study of such systems, as well as other examples such as 
chemical pattern-forming reaction-diffusion systems, has also revealed that 
an accurate description of their dynamics may require inclusion of the effects 
of spatially-inhomogeneous distributions of molecules. The usual framework 
to analyse such situations is that of stochastic reaction-diffusion models. 

A popular approach for studying stochastic reaction-diffusion models in¬ 
volves decomposing the domain occupied by the system into small compart¬ 
ments or voxels, and counting the number of molecules of each species in 
each voxel [331 U- Chemical reactions between species are treated locally 
(i.e. within each voxel), their rates being determined by the local abundance 
of each constituent species. Molecular transport between compartments is 
usually modelled by simple diffusion. 

This compartmental approach poses a problem: numerical (Monte Carlo) 
techniques such as the Gillespie stochastic simulation algorithm f SSA) [T8lfT^ 
or the r—leap method [201 E] are inefficient since they scale poorly with 
the number of reaction channels, which is proportional to the number of 
compartments. 

In order to overcome issues regarding the efficiency of direct simulation 
methods, several strategies have been proposed. For example, the Next Re¬ 
action Method (NRM) proposed by Gibson & Bruck [T7[ is an exact method 
in the same sense as the SSA, i.e. the sample paths generated are exact 
realisations of the solution of the corresponding Master Equation. However, 
unlike the SSA, the computing time grows logarithmically with the number 
of reaction channels. This is accomplished by (i) recycling the random num¬ 
bers generated in previous time steps, so at each time step only one new 
random number must be generated, and (ii) organising the reaction channels 
in a queue (more specihcally, a tree) with events ordered in ascending order 
of waiting time. In this way, it is possible to determine which channel will hre 
next when by looking at the root of the tree. The Next Sub-Volume Method 
(NSVM) is a modihcation of the NRM which accounts for reaction-diffusion 
systems ra. It deals with the stiffness problem associated to the fact that 
the transition rates associated with diffusion are proportional to h~'^, where 
h is the compartment length. 
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An alternative set of methods are hybrid methods. Often a stochastic 
description is only needed in a certain region of the domain. Elsewhere 
the number of components is large enough that a mean-held description is 
reasonable. A paradigmatic example of this situation is front propagation in 
reaction-diffusion systems (HEIE]- In systems such as the stochastic Fisher- 
Kolmogorov model, the number of particles ahead of the front is small and, 
therefore, huctuations need to be taken into account. By contrast, behind the 
front, the number of particles huctuates about the carrying capacity of the 
model. In these conditions, simulating the system using the compartmental 
approach and a direct (Gillespie) simulation method is inefficient and it is 
natural to propose a hybrid approach where the region behind the front 
is modelled as a Fisher-Kolmogorov PDF, the region ahead of the front is 
treated as a stochastic process, and appropriate matching conditions are 
applied in the intermediate region. 

This idea, and variations thereupon, have been implemented for several 
different methods and several systems. Its hrst incarnations consisted of 
hybrid models for pure diffusion [T^ flTl 115] . In [15] uni-molecular chemical 
reactions such as chemoabsorption are also considered. In all these methods 
the boundary or overlapping region between the two regimes is considered 
hxed. 

A further step forward towards algorithms that can cope with more com¬ 
plex situations has recently been proposed by Hellander et ah |2T] . Based on 
previous work on a method to simulate stochastic react ion-diffusion systems 
on unstructured domains nn, they have proposed a method where, at each 
voxel, the different species are divided into two classes, namely, mesoscopic 
(i.e. well-mixed within the voxel and with dynamics determined by the cor¬ 
responding Master Equation) and microscopic. Concerning the latter, they 
are assumed to be subject to off-lattice reaction-diffusion dynamic, modelled 
in terms of the Green’s function method for the corresponding Smoluchowski 
equation proposed by Van Zon and ten Wolde [371ES] . The coupling between 
mesoscopic and microscopic degrees of freedom within each voxel is accom¬ 
plished via a splitting scheme. Similar hybrid algorithms are discussed in 

PElEHEHlIISllaollSIlESlEHI. 

Moro [5] proposed a similar method for simulating stochastic reaction- 
diffusion systems with propagating fronts in which a macroscopic PDE is 
coupled with a mesoscopic Master Equation. The two descriptions hold in 
different sub-domains and are coupled across a moving boundary, using a 
method which balances the fluxes between the sub-domains on average only. 
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The hybrid algorithm presented in this paper also couples a mesoscopic 
description of a stochastic reaction-diffusion system, modelled by an on- 
lattice, react ion-diffusion Master Equation (RDME), with a reaction-diffusion 
system, which is obtained from the mean held equations associated with the 
stochastic model. We remark that whilst our mean-held reaction-dihusion 
system converges to a reaction-dihusion PDE in the continuum limit, the 
RDME does not converge to the associated oh-lattice (Doi or Schmolokowski) 
reaction-dihusion models |2Tl[22] in spatial dimensions larger than 2. In light 
of this, our method should be viewed as the hybridisation of the RDME with 
its associated on-lattice, mean-held reaction -dihusion limit. 

Our method extends the work of [S] to systems without propagating 
fronts. Furthermore, the interface condition used here preserves the total 
amount of particles at all times in every single simulation. The interface is 
also chosen such that the mean held description is still valid at the interface 
region, and each interface is only one compartment in size. This way, errors in 
the hux between the two domains are negligible. We allow multiple interfaces 
so the stochastic and deterministic regions need not be connected, and inter¬ 
faces may be dynamic in space and time. For multi-species models, diherent 
species may exist in diherent stochastic and deterministic sub-domains. We 
test our algorithm on two classical systems, the Fisher-Kolmogorov equation 
and a spatial Lotka-Volterra system. The former serves as the simplest ex¬ 
ample of a system with a propagating front, and in our hybrid model a single 
interface separates the stochastic region at the wave front from the deter¬ 
ministic region behind the wave front. The spatially resolved Lotka-Volterra 
system serves as a test case for a multi-species reaction-dihusion system, 
where the hybrid model can, in general, have multiple species-specihc inter¬ 
faces. 

The paper is organised as follows: In section we introduce notations and 
conventions for stochastic and deterministic reaction-dihusion systems. The 
methodology for a generic hybrid model with an arbitrary number of species is 
presented in section and our algorithm is applied to the Fisher-Kolmogorov 
equation in section]^ Finally, the spatially resolved Lotka-Volterra model is 
investigated in section 
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2. Stochastic and Deterministic Reaction-Diffusion Systems 

2.1. Stochastic Reaction-Diffusion Systems 

We consider a system of Smax species defined on a regnlar lattice. Each 
species comprises individnals that can migrate to neighbonring lattice sites, 
or react locally with entities of the same or other species. For simplicity, 
we describe the system in one spatial dimension, where the lattice index 
is k = 1,..., kmaxi hnt the generalisation to higher dimensions is straight¬ 
forward. We let Ns{k,t), s = 1,... ,Smax denote the nnmber of individnals 
of species s in box k at time t, and let h denote the lattice constant, so 
that L = hkmax is the domain size. If Ns{k,t) is Markovian, then the time 
evolntion of Ns{k,t) is governed by a master eqnation of the form 

= ( 1 ) 

N 

Here, N denotes a generic state specified by the nnmber of individuals Ns{k, t) 
of all species s in any compartment fc, and likewise iV, so the sum is under¬ 
stood to be over all states defined by N. The probability that the system is 
in state N at time t is denoted P{N,t) and the transition rate for a change 
from state N to state N is In what follows, we suppress the time de¬ 

pendence of Ns{k,t) whenever it is clear that Ns{k) depends on the current 
time t. We will consider two types of transition rates one describing a 

random walk and the second local reactions: 

TN,(k)-l,Nsil)-i-l\Nsik),N,il) = -J^^s{k), I = k±l 

'^Ni(k)+pi^r,-,Nsraaxik)+Psmax,r\l^l{k),..;Nsmaxik) ~ {N,{k),...,N,_{k)). (2) 

In k = 2,..., kmax — 1; wlicreas the specification of the transition rates 
in the boundary boxes k = l,kmax depends on the boundary conditions. 
Furthermore, r = 1,..., M denotes the index of a particular reaction, and 
Ps^r specifies how reaction r changes the number of individuals of species 
s. Thus, a step in the random walk changes the spatial distribution of a 
particular species, whereas reactions act locally in a particular box k. 

We define shift operators 

EtJ(N,(k),...)-f(N,{k)±l,...), (3) 
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where the dots indicate that the function can depend on Nsi{k') for s' ^ s, 
k' 7 ^ k, but operator affects only Ns{k). Substituting ([^ in Q, we can 
write the master equation as 


dP 

dt 



2.2. Mean-Field Limit 

The mean quantities Ns{k) = Ns{k)P{N) evolve in time according 



r 


(5) 



The equations for fc = 1, kmax depend on the boundary conditions. If the 
reaction terms are non-linear, then these equations are not closed but de¬ 
pend on higher moments of Ns{k). However, if the number of particles of 
each species is large, Ns{k) oc H, where H 3> 1, we can perform a system 
size expansion in inverse powers of H [35]. The leading order term in the 
expansion closes the time evolution equations for the means so that 



N,{k + 1) + Ns{k - 1) - 2N,{k) )+RA iVi(fc),..., , 


( 6 ) 


By expanding Ps,rRr {{Ni{k ),..., to lowest order in jt, we ob¬ 

tain the total reaction rate Rg for species s. In the limit h —)■ 0 we can 
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obtain from ([^ continuum react ion-diffusion equations for the particle den¬ 
sities ns{x, t): 


= D,V^n,{x, t) + r,(ni,..., n,_J. (7) 

In ([^ we have abused notation, identihed ns{x,t) = ns{k,t) = for 

kh — ^ < X < kh + ^, and introduced a focal reaction term r^, which is 
obtained from Rg by 


Rg ( N^{k), ..., ) = hvg 


( iVi(fc) 




( 8 ) 


Note that in higher spatial dimensions {d > 2), the h in the denominator 
needs to be replaced by Furthermore, we note that to solve Q numeri¬ 
cally, we must discretise the PDF, potentially with a different discretization 
than that used for the stochastic compartment model. This might be de¬ 
sirable as for a PDF, typically we want the lattice to be as small as com¬ 
putationally reasonable to avoid discretization errors; by contrast for the 
stochastic model we must be careful when simulating nonlinear reactions if 
the lattice becomes similar in size to the reaction radii. Hence, if we do not 
alter the compartmental model, we cannot decrease the compartment size 
arbitrarily, see also the discussion in [2l]. For the remainder of this paper we 
use the mean held equation ([^, and do not consider the continuum limit Q 
any further. We emphasise wheather a variable is part of the deterministic 
or stochastic regime by using capital letters for stochastic variables, and 
lower case letters for deterministic variables, so (|^ is rewritten in terms 
of the Uk as 

^ {ng{k + 1) + ng{k - 1) - 2n,(/c)) + r, (^ui(/c),..., , 

k 2,. .., kjYiax !• (9) 

2.3. Method of Solution of the Stochastic Model 

We solve the stochastic model using the Gillespie algorithm [THl [I9] , ex¬ 
ploiting the fact that the time to the next event is distributed exponentially. 
Hence, if denotes one of the non-zero transition rates of the model, and 

















Ti is a random number uniformly distributed in [0,1], then 


T = — 



( 10 ) 


gives the time to the next event. Furthermore, if r 2 is a second, independent 
random number uniformly distributed in [0,1], we can calculate which event 
I happens by imposing the condition 



( 11 ) 


The advantage of using the Gillespie algorithm is that it simulates exactly the 
dehned Markov process. The hybrid methods discussed in this paper do not 
depend on the details of how the stochastic process is simulated, and faster 
algorithms such as the r leaping algorithm 120118] can also be used. Since 
these algorithms are not exact, in this paper we prefer to keep the stochastic 
simulation as detailed and accurate as possible and use the original Gillespie 
algorithm, and improve computational performance by switching to a mean 
held description when appropriate. 

3. General Algorithm for Stochastic/Deterministic Reaction-Diffusion 

Hybrid Model 

Our algorithm involves decomposing the spatial domain into two regions 
for each species. In one region the species of interest is modelled in a stochas¬ 
tic way; in the other the mean-held limit is used. In regions where one 
species is modelled stochastically and the other deterministically, all interac¬ 
tion terms involving the two species are modelled stochastically. The inter¬ 
face condition describes how the two domains are coupled together, and how 
the interface moves. The following steps are performed by the hybrid model 
repeatedly: 

Hybrid Algorithm 

1. Generate time to next stochastic event via equation ( [l0| ) 

2. Simulate which stochastic event happens via condition (0 

3. Iterate hnite diherence scheme to new time via equation ([^ 

4. Galculate interface condition 
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5. Return to Step 1 until specified end time 


For notational simplicity, in the main part of this section, we focus on a 
one-species system with a single interface. We explain in Appendix A how 
this algorithm can be modihed to account for multiple species and multiple 
dynamic interfaces. 


3.1. Hybrid Algorithm with a Single Interface 

Let [0, L] be our modelling domain, discretized as before into k^ax com¬ 
partments of size h such that k = 1,..., fc/ — 1 is the mean held domain, k = 
kj + 1,..., kmax is the stochastic domain, and k = kj labels the interface com¬ 
partment. We use equation ([^ to solve for the variables ns(l),... ,ns(fc/_i) 
in the deterministic regime, whereas the evolution of the stochastic variables 
Ns{ki+i),... ,Ns{kmax) is determined by simulations of the master equation 
([^ with transition rates ([^. The equations used to determine the variables 
at kj will be discussed below. 


Mean Field 

Interface 

Stochastic 

Deterministic Reactions 

Stochastic Reactions 

Stochastic Reactions 

Deterministic Flux 


Stochastic Flux 


Figure 1: The domain is decomposed into a deterministic region where the 
system is described by the mean held equations, and a stochastic region in 
which it is described by the stochastic equations. At the interface, which is 
a single compartment, between these domains the hux into the mean held 
domain is deterministic (equation ([T^), whereas the reactions and hux into 
the stochastic domain are calculated in a stochastic way (equations and 

(P))- 


3.1.1. Fluxes and Reactions at the Interface 

We now explain how the stochastic and deterministic regimes are coupled 
at the interface. We identify at all times 


ns{kj) 


Nsjki) 

h 


( 12 ) 


10 














to emphasise that the interface compartment will exhibit both deterministic 
and stochastic behaviour. Three processes contribute to changes in particle 
numbers in compartment kj: fluxes into and from compartment {kj — 1), 
which is part of the mean held domain, huxes into and from compartment 
(/c/ + 1), which is part of the stochastic domain, and local reactions (see 
Figure]^. Hence, we model the hux between compartments kj and {kj — 1) 
deterministically, and the hux between compartments kj and {kj + 1) in a 
stochastic manner. If r denotes the current Gillespie time step, we calculate 

Usiki, t + r) = ns{ki, t) + (n^(/c7, t) - ns{ki - 1, t)). (13) 

The hux between boxes kj and (fc/ + 1) is accounted for by the transition 
rates 


TNs{ki)-l,Ns{ki+l)+l\N(ki),N{ki+l) — J^Ns{ki), 

TN4ki)+l,Nsiki+l)-l\Nsiki),Nsiki + l) = J^^s{kl + 1). (14) 


We also have 


TN,(ki}±l,Nsiki-l)^l\Nsiki),Nsiki-l) — 0 , 


(15) 


as the corresponding hux is already accounted for by equation (13). Finally, 


we specify the local reactions in a stochastic way via transition rates 


% 


Nl{kl)+pi,r,---,Nsmax{^l)+Psmax,r\Nl{kl),...,Ns, 


( 16 ) 

In [Appendix A.3 we discuss an alternative formulation for which local reac¬ 
tions are calculated in the mean held framework. 

We remark that due to the interfacial coupling, the mean held solution 


acquires some stochasticity. Indeed, formally, (13) appears to correspond to 
a Neumann-like boundary condition at the interface. However, since Ugiki) is 
subject to both stochastic reactions and a stochastic hux into the stochastic 
domain, the interface condition appears as a stochastic source for the mean 
held model at the interface. 

The mean behaviour of the full hybrid model is obtained by solving the 
mean held equations in the whole domain. By construction, these are ob¬ 
tained from the mean of the stochastic model. We now have to convince 
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ourselves that the mean behaviour at the interface gives the mean held limit 


of the full stochastic model. Note that (13) appears to diverge in the limit 
h —)■ 0. This is because the deterministic contribution due to diffusion in¬ 
cludes only the hux between the interface and the deterministic domain. 
To obtain the full discrete Laplacian, we must add the mean hux to the 
right-hand side. This is easily obtained if we calculate the mean held limit 
associated with the transition rates (14), giving ^ {ns{ki -|- l,t) — ns{ki,t)), 
and this contribution is simply added to (13). 


Fractional particles. Conventionally, the state of the stochastic model is de- 
hned by the numbers of individuals Ns{k) in each compartment k, and these 
are non-negative integers, whereas the densities ns{k) are real valued. At the 
interface the hux into and from the mean held domain is given by equation 


(13) which alters ns{ki), and thus Ns{kj), by a real valued number. However, 
local reactions, as well as the hux into and from the stochastic domain, are 


described by transition rates (14) and (16), which ehect integer changes in 
Ns{ki). Consequently, it is not a priori clear whether the stochastic com¬ 
ponents of the reactions and huxes at the interface are well dehned. First, 
we note that, by dehnition, the interface is such that the particle number 
there is sufficiently large, Ns{ki) S> 1, so that the mean held description is 
accurate, and hence agrees closely with the corresponding stochastic model 
(which has integer valued particle numbers). Adding a real part between 
zero and one to a large integer number will not signihcantly alter the transi¬ 
tion rates, so the stochastic model with real valued Ns{ki) will agree closely 
with both the stochastic model with integer valued Ng{kj) and the mean 
held model. Formally, the state space of the stochastic model at the inter¬ 
face is thus still in one-to-one correspondence with the integers, which are 
shifted by the fractional part of Ns{k). Hence, the stochastic part of the hy¬ 
brid model is well-dehned. There is still a numerical mismatch between the 
results obtained from a Gillespie algorithm with fractional or with integer 
valued numbers, but as long as Ns{ki) S> 1, we found this mismatch to be 
negligible. 


3.2. Moving Interface Condition 

The condition used to locate the interface is not dictated by a rigorous 
mathematical requirement: it represents a compromise between performance 
and accuracy. We view the mean held equations as an approximation to the 
stochastic model that neglects huctuations. Hence, the larger the mean held 
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domain, the fewer stochastic fluctuations are taken into account. However, 
we then typically increase the performance of the simulation. 

We determine whether a compartment belongs to the stochastic or mean 
held domain by comparing the number of particles in that compartment with 
a threshold number, 0^, such that if Ns{k) < 0^, then box k is part of the 
stochastic domain, and otherwise part of the mean held domairj^ This is 
justihed since huctuations typically scale with the square root of the number 
of particles in a box, ^^/N^(k). As the number of individuals may change over 
in time, the position of the interface may also evolve in space and time. 

We also implement a minimum domain size condition as a simple check 
that no connected component of the mean held domain is allowed to become 
too small. Imagine, for instance, that the mean held domain is enclosed by 
two disconnected components of the stochastic domain. If the mean held 
domain comprises only a few compartments in the discretisation, it might be 
computationally more efficient to remove the mean held domain and absorb it 
into the stochastic domain, as then we do not have to calculate the interface 
condition. In the simulations performed in this paper a minimum domain size 
of 5 compartments for the mean held model was used. In higher dimensions, 
we anticipate that a cube with a length of 5 compartments would work well. 
We stress, however, that the choice of threshold conditions is model specihc, 
and hence the minimal domain size requirement should be chosen on a case- 
by-case basis. 

After the position of the interface is updated, we must check that all 
particle numbers in the stochastic domain are integer values, paying partic¬ 
ular attention to compartments that were previously part of the mean held 
domain. By mass conservation, this will result in a renormalisation of the 
density functions in the mean held domain: this procedure is discussed in 
section 13.31 

The following steps are used to adjust the position of the interface: 

1. Calculate threshold condition to locate interface; 

2. Calculate minimum domain size condition; 

3. Renormalise particles and densities. 


^If this condition splits the whole modelling domain such that the stochastic domain 
consists of multiple, disconnected components we either need to introduce multiple inter¬ 
faces, seelAppendix A. 41 or connect the disconnected components. 


13 





3.3. Renormalision of Particle Distribution 

When the interface moves such that a compartment previously treated 
deterministically and, hence, described by real valued densities ns{k), enters 
the stochastic domain, we need to ensure that the number of particles be¬ 
comes integer valued. By mass conservation, we cannot simply remove the 
fractional part of Ns{k) = hns{k)\ instead we rescale the densities outside 
the stochastic domain. 

Let us focus on a single box k (and assume that the single interface has 
moved by only one compartment, so that Ng{k) is now non-integer valued 
but part of the stochastic domain). We interpret the fractional part 

Ps = Ng{k) mod (1), (17) 


as the probability that an additional particle is in this box. We draw a 
uniform random number r G [0,1]. If r < ps then we place the particle 
in box k] otherwise it is placed in the deterministic domain. To preserve 
particle numbers, we renormalise the density function in the deterministic 
regime and reset Ns{k) so that if r < ps then 


Ng{k) ^Ng{k) + {l-pg), 




ng{l), l = l,...,k-l, 


(18) 


and otherwise 

Ng{k)^Ng{k)-pg, 

n,(/)^(l+ [ n,(/), l = l,...,k-l. (19) 

A similar rescaling procedure was used in [15], where a reaction-diffusion 
PDE was coupled with a Brownian dynamics model. There, it was found 
that particles crossing the interface twice can cause increased variance. Our 
algorithm does not lead to an observable increase in variance as, by con¬ 
struction, the interface is chosen so that the mean field and interface domain 
always contain a large number of particles. 
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4. Stochastic Fisher-Kolmogorov Equation 

In this section, we apply onr algorithm to the Fisher-Kolmogorov eqna- 
tion, 


On d‘^n , _ n. 

= + And - jj) 


( 20 ) 


Here, D is the diffusion coefficient, A the growth rate and lA the carrying 
capacity. 


A stochastic, lattice-based version of equation (20) was studied in m 
and is defined by a master equation with the transition rates 

T/V(fc)-l,7V(fc±l)+l|Ar(fc),Ar(fc±l) = 

TN{k)+l\N{k) = XN{k), 

TN{k)-l\N{k) = -^N{k){N{k) - 1), k = 1,. . . ,kmax, (21) 

with Dirichlet boundary conditions Ni = = 0. The evolution of the 

means N{k) is given by 


dN{k) D 


dt 


= ^ + 1) - SiV)*;) + N{k - 1)) + \N{k) fl - (22) 


Defining n(x,t) = with x = hk, the mean field equation for n reduces 


to (20) if D 3> 1, as the van Kampen approximation at leading order implies 

'-' - - 2 

the moment reduction N{k,t)‘^ N{k,t) . 

We will now use the hybrid algorithm to simulate travelling wave solutions 
as we vary several model parameters. As noted in HI M for the stochastic 
Fisher-Kolmogorov equation with the conventions used in the present paper 
and in for alternative formulations, stochastic effects can produce 

wave speeds Cstoch which deviate from the deterministic Fisher-Kolmogorov 
equation (20), cpde = 2\/DX. The three other parameters control the var¬ 
ious limits: 0 —?• cx) yields the stochastic model from the hybrid model, 
D —>■ cx) yields the mean field model from the stochastic model, and h —0 
yields the PDF from the mean field model. We will study the effect of 
variation in these three parameters on travelling wave speeds in the next 
subsection. 
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4.1. Fisher-Kolmogorov Travelling Waves 

We now study travelling wave solutions, fixing D = 1, A = 1. We ensure 
that in all simulations the travelling wave is sufficiently far from the bound¬ 
aries, so that effects associated with the hnite domain size are negligible. 
As initial conditions we approximate the travelling wave solution of the PDE 


(20)), so that we can focus on wave propagation, rather than wave formation. 


Figure [^compares travelling wave solutions generated from the stochastic 
model, the hybrid model (with 0 = 25) and the PDE. In each plot we 
present a single realisation and the mean of 256 realisations of the stochastic 
(Figure 2(a)[|(c)[(e) ) and hybrid models (Figure 2(b)[(d)[(f) ) together with 
the numerical solution of the corresponding PDE. We £x D = A = 1, kmax = 
20, h = 1, and allow to vary D. We note that the travelling wave speeds for 
the stochastic and hybrid models are slower than those of the PDE, and the 
speed increases with D. Furthermore, the relative noise, i.e. the fluctuation 
of a single stochastic or hybrid realisation about the mean, decreases as D 
increases. Finally, the wave front of the PDE appears to be steeper compared 
to the wave front of the mean of 256 realisations both in the stochastic and 
hybrid models compared to the PDE, and the steepness increases with D. 
This is explained as the different realizations of the stochastic model can 
have different speeds, hence the average broadens the wave front. 

We now compare the stochastic to the hybrid model. For D = 10, so 
D < 0, neither single realisations nor the mean of the stochastic model 


(Figure 2(a)) differ signihcantly from the hybrid model (Figure 2(b)), as 


the threshold of the hybrid model is considerably larger than the carrying 
capacity, so almost certainly the entire domain of the hybrid model will be 
stochastic. When D = 25, so D = 0, the stochastic model (Figure [2(c)[ ) and 
the hybrid model (Figure 2(d)) differ signihcantly away from the wave front. 
The stochastic model is much noisier, but the noise in the hybrid model is 
non-zero as noise from the stochastic domain can diffuse into the mean held 
domain, raising the particle number above the threshold value. Note that 
huctuations can also reach beyond the threshold in the hybrid model due to 
the minimum domain size requirement for the deterministic domain. When 
D = 50, so D > 0, the noise a way from the wave front associated with the 
stochastic model (Figure [2(e)[ ) is absent in the hybrid model (Figure [2(f)| ). 
Nevertheless, the wave fronts of the means appear similar. Hence, for the 
parameter values used, the hybrid model represents a good approximation 
to the stochastic model, producing travelling waves with the same speed. If 
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the carrying capacity is larger than the threshold 0 then fluctuations around 
the carrying capacity behind the front are suppressed, without affecting the 
wave speed. 

Figure l^shows the dependence of the wave speed on the lattice constant. 
We compare the stochastic model against the mean held model, i.e. the 
finite difference discretization of the PDF with the same lattice constant, 
and the hybrid model. Densities are fixed by adjusting Q = 80* h. Likewise, 
0 = 16 * h, but we also compare to the hybrid model with a fixed threshold 
0 = 10. The wave speed c is calculated by observing that the change of 
the total number of particles in time, averaged over all simulations, Wot = 
YllZT should be proportional to c. We approximate the wave speed 

by comparing Wot after fixed time intervals At = 5, obtaining Ntot{t + At) — 
Ntotit) = The finite difference model converges, as expected, to c = 2 

as h —)■ 0, even though every finite lattice spacing will still result in some 
visible dispersion after a long time. The stochastic model is significantly 
lower the wavespeed while the wave speed of the hybrid model, with 0 = 10, 
converges to that for the stochastic model. However, if 0 = lOh, then the 
hybrid model does not appear to converge to the stochastic model. 

Figure 4(a) shows how the wave speed varies with the carrying capacity D 
for the stochastic and hybrid models. Here, h = D = A = 1 in all cases. The 
wave speeds for the stochastic model are identical to those obtained in | 1 ] , and 
for large D they approach the wave speed of the mean held theory as expected 
by general theory [26l[35], this value is slightly above c = 2 because h is finite. 
When 0 = 100, the wave speeds for the hybrid and stochastic models are 
indistinguishable. We remark that a naive expectation that the hybrid model 
should be intermediate between the PDF and the stochastic model does 
not imply that the wave speed of the hybrid model should be intermediate 
between their wave speeds. However, a correct expectation is that as 0 
increases, the agreement between the hybrid and stochastic models, including 
their wave speeds, increases. Figure 4(b)| shows this clearly. Here, we have 
fixed D = 100 and varied 0. The wave speed monotonically approaches that 
of the stochastic model (not shown, but it is formally obtained as 0 —)■ oo), 
and for 0 ~ 25 the wave speed of the hybrid model is almost identical to 
that of the stochastic model. Due to the observation in those numerical 
simulations that the wave speed seems to plateau for high and low values of 
0 , we have fitted the function c = Oi + a 2 Erf {a^*Log{Q) + 04 ) to the values 
obtained from the simulations, and obtained a very good fit for the values 
ai = 1.5, 02 = 0.13, Os = 1.4, 04 = —3.0. 
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(c) Stochastic model, O = 25 



200 


(d) Hybrid model, 0 = 25, fl = 25 




200 


(f) Hybrid model, 0 = 25, O = 50 

Figure 2: Series of plots comparing the travelling waves profile generated by 
the stochastic model in the column to the left and the hybrid model with a 
threshold of 0 = 25 on the right (simulation time t = 60). As the carrying 
capacity increase, 12 = 10, 25,100. Each plot shows a single realisation 
as well as the mean of 256 realisations of the stochastic or hybrid model, 
respectively. Results for the corresponding PDE (20) are also shown. The 


other parameter values are = A = 1, kmax = 20, h = 1. 
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Lattice Spacing h 

Figure 3: Series of curves showing how the wave speed of the stochastic 
model, the mean held model (i.e. the hnite difference discretisation of the 
PDF) and hybrid models for thresholds of 0 = 10 and IQ * h changes as 
the lattice spacing varies. We note that the wave speed of the hybrid model, 
with a hxed threshold, converges to that of the stochastic model, whereas the 
hybrid model where the threshold is adjusted with the lattice spacing does 
not. The other parameters are r2 = 80*/i, Z1 = A = 1, and all stochastic and 
hybrid results are obtained from averaging 1024 different simulations. 
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Figure 4: The average wave speed dependence on the carrying capacity is 
shown for the stochastic model, as well as the hybrid model for thresholds 
of 0 = 10, 25 and 100 in (a) The dashed-dotted graph of the hybrid model 


with a threshold of 0 = 100 coincides with the solid line representing the 
stochastic model. The wave speed of the PDF is c = 2, but the discrete mean 
held equations can slightly deviate from this value in an ^-independent way. 


(b) shows, for a hxed carrying capacity of = 100, explicitly the threshold 


dependence for several values of 0 , and the solid line is a ht via the function 
c = ai + a 2 Erf{a 3 * Log{Q) + 04 ), where Erf is the error function and we 
obtained Oi = 1.5,02 = 0.13,03 = 1 . 4,04 = —3.0. The other parameters 
are h = D = X = 1, and all stochastic and hybrid results are obtained from 
averaging 256 different simulations. 
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5. Spatial Stochastic Lotka-Volterra System 

We now investigate a predator-prey system predator M and prey N. In 
the stochastic model, each species can jump to neighbouring lattice sites, 
the prey reproduce at rate a, predators die at rate c and consume prey and 
reproduce at rate b. The transition rates are given by 


TN{k)-l,N{k±l)+l\N{k),N{k±l) 

TM{k)-l,M{k±l)+l\M{k),M{k±l) 

TN{k)+l\N{k) 

TN{k)-l,M{k)+l\N{k),M{k) 

TM{k)-l\M{k) 


aN{k), 

hN{k)M{k), 

cM{k). 


(23) 


For simplicity, we choose the interaction reaction such that each time a prey 
is eaten by a predator, a single new predator is born. The mean held and 
continuum limit corresponds to the classical spatial Lotka-Volterra equations: 


dn 

di 

dm 


Dn 

Dm 


d‘^n 

d'^m 

dx'^ 


+ an — bnm, 
-|- bnm — cm. 


(24) 


Here, n = n(a;, t) and m = m(a;, t) are prey and predator densities related to 
N{k) and M{k) repectively in the same way as before. 

As before, we use a hnite difference approximation to solve equation (24), 
discretising in space using the same lattice as for the stochastic model. For 
the time integration we use the Runge-Kutta method. All plots are nor¬ 
malised so that the number of predators or prey in a box is shown, rather 
than the corresponding density. 

There are now four subdomains to consider, depending on whether each 
of the predator and prey evolve deterministically or stochastically. For no- 
tational simplicity we identify N{k) = hn{k) and M{k) = hm{k) where 
appropriate. We will now explain explicitly which transition rates dehne the 
stochastic model, and which PDFs correspond to the mean held model solved 
in the respective subdomain. The interfaces between the subdomains are as 


given in section 3.1.1 
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Deterministic Predator and Deterministic Prey 

In this region, we solve Equations (24), and set the transition rates (23) 
equal to zero, so no stochastic reactions occur. 


(26) 


Deterministic Predator and Stochastic Prey 
The deterministic equation are: 

dn dm d^m 

^ = 0 . 

and the transition rates are the following 

'TN(k)-l,N{k±l)+l\N(k},N(k±l) = {k) ,TM{k)-l,M{k±l)+l\M{k),M{k±l) = 0, 

Tv(fc)+l|7V(fc) = 0,N{k),TN{k)-l,M{k)+l\N(k),M(k) = bN{k)M{k),TM(k)-l\M(k) = 0. 

(26) 

Stochastic Predator and Deterministic Prey 

Deterministic part equations are described below. 


dn d'^n dm 


(27) 


and the transition rates are the following 


D 


TN{k)-l,N{k±l)+l\N{k),N{k±l) — 0, 7M(fc)-l,M(fc±l)+l|Ar(fc),M(fc±l) — 

TN(k)+l\N(k) = 0, 7iV(A:)-l,M(fc)+l|A7(fc),M(fc) = bN {k)M (k) , TM{k)-l\M{k) = cM{k). 

(28) 

Stochastic Predator and Stochastic Prey 


Here both species are fully stochastic and we use transition rates (23). 


We will now consider two scenarios appearing in the spatial Lotka-Volterra 
system, and compare the hybrid model to the stochastic and deterministic 
models. Both scenarios correspond to solutions of the PDE which oscillate 
in space and time, but in one case the oscillations bring the total number of 
prey so close to zero that extinction is possible in the stochastic model. 
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5.1. Oscillatory Behaviour without Observable Extinction 

The domain of length L = 20 is divided into kmax = 101 boxes, so h = 0.2. 
Initially it contains a spatially homogeneous distribution of prey and preda¬ 
tors so that N{k,t = 0) = 50, M{k,t = 0) = 5. The model parameters are 
hxed so that Dj^f = Dm = 1, a = 2, 6 = 0.1, c = 3. With Neumann boundary 
conditions, equations (24) remain spatially homogeneous at all times, and 
both populations oscillate in time. Typical results are presented in Figure 



Figure 5: The solution of the Lotka-Volterra Eq.(24) with parameters 
a = 2,b = 0.1, c = 3 for spatially homogeneous initial conditions N{k,t = 
0) = 50, M{k,t = 0) = 5. Shown is the time evolution of the number of prey 
and predators in any given box k in the discretisation, to allow for better 
comparison with the stochastic model. As the diffusion terms do not con¬ 
tribute in the spatially homogeneous case, this solution is identical to the 
solution of the Lotka-Volterra ODEs. 


1^ show that the peak in prey numbers is followed by a peak in the number 
of predators. For this choice of parameter values the minimum number of 
individuals of either species is always sufficiently large that extinction in the 
stochastic, spatial model is almost impossible. Corresponding results for the 
stochastic and hybrid models for two choices of the threshold values (0 = 10 
and 0 = 25) are shown in Figure]^ The column on the left shows the spatial 
profile of the number of predators and prey in a given box at time t = 2.5, 
both for a single realisation and the mean of 256 different realisations. Com¬ 
paring either predator and prey numbers, we note that the means for the 
stochastic and hybrid models appear similar for both values of 0 (see Fig- 
6(a)|p(c)|[(e) ), and are fairly homogeneous, whereas single realisations of 


ure 
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either model differ markedly. As expected, the profile of predator and prey 
numbers in the stochastic model is noisy throughout the domain, whereas 
noise is suppressed in the hybrid model when the population numbers ex¬ 
ceed the threshold. We observe the predator numbers in Figure |6(c) are 
above the threshold and hence smoothly distributed in space, but they are 
not homogeneous, in constrast to the profile of the mean. 

The column to the right in Figurej^shows the time evolution of the spatial 
mean of the number of prey, {N)^^ = T~Ylk=T^i^)- 
means for the stochastic model (Figure 6 (b")| and the hybrid model are similar 
for either threshold (Figures 6 (d)[[(f) ). In all cases, we observe oscillatory 
behaviour around 30 = |, the steady state of the corresponding PDF Eq. 
(24), with a decreasing amplitude. This contrasts with the dynamics of 
the PDE (see Figure]^ where the amplitude of oscillations was constant in 
time. This is because stochastic fluctuations may cause oscillations to fall 
out of synchrony and hence oscillations average out. Hence, the damping 
effect is stronger when plotting the mean of all realisations, rather than a 
single realisation, where the amplitude can fluctuate in time. However, this 
also implies that the PDE is not a good approximation to the mean over 
different realisations. We finally remark that while the plots of the spatial 
mean number of predators look different for the stochastic and hybrid models 
with the two thresholds, this is not significant as different realisations of the 
same model (not shown here) also look different. To properly compare the 
stochastic and hybrid models, we need quantitative measures (Table [^. We 



Stochastic 

Hybrid 0 = 25 

Hybrid 0 = 10 

Mean Field 

(N),, 

30.50 ±0.05 

30.32 ±0.05 

30.30 ±0.16 

30.47 


20.35 ±0.04 

20.23 ± 0.04 

20.22 ±0.13 

20.20 

d wl, - ml. 

16.0 ±0.3 

16.0 ±0.3 

15.4 ±0.4 


- (M)l, 

12.8 ±0.3 

12.8 ±0.3 

12.3 ±0.4 



Table 1: Mean of prey population over space, time and 256 different real¬ 
isations, Eq. (29), as well as the corresponding number of predators. The 

standard deviation sj — {N)j^ ^ = sj ^ X]r=i of prey numbers, 
and likewise for predator numbers, are also shown, as are results for the cor¬ 
responding mean field model, where there are no different realisations, and, 
hence, no standard deviation. 
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calculare the spatial and temporal average of the mean nnmber of prey across 
256 different realisations, 


256 K 




‘25Qkffiaxtmax 


EE 

r=l k=l 


Nr{k, t)dt, 


(29) 


over the simnlation time of t^ax = 50, and likewise the average of the mean 
nnmber of predators. In this way we obtain a single nnmber with large sta¬ 
tistical signihcance, which is easier to compare as at individnal time points, 
the oscillations in different realizations can be ont of synchrony. We observe 
good agreement between the models. The valnes reported in Table are 
in good agreement with the steady state valnes of the corresponding ODE 
model, bnt slightly different, as the temporal oscillations are not necessarily 
symmetric with respect to the steady state valne. We then measnre the stan¬ 
dard deviation of the spatio-temporal averages, which is — (iV)^ ^ for 

prey nnmbers, and likewise for predator nnmbers. Here, we note that the 
stochastic model is in close agreement with the hybrid model with 0 = 25; 
this agreement is less striking when 0 = 10. We conclnde that the hybrid 
model can reprodnce stochastic measnres of the stochastic model, snch as 
the standard deviation, and the agreement is better for larger valnes of the 
threshold 0. 


5.2. Extinction and Blow-up 

We now choose a spatially homogeneous population of prey, N{k, t = 0) = 
50, A; = 1,..., kmax, a number of predators present only on the left side of the 
domain, M{k) = 100, (fc = 1,..., 9), M(10) = 98, M(ll) = 50, M(12) = 2, 
and M{k) = 0, (fc = 13,... 101), Neumann boundary conditions and param¬ 
eters Dtv = Dm = 1, a = 1, 5 = 0.1, c = 2, L = 20, kmax = 101. This scenario 
is typical example of invasion of a predator into a population of prey, leading 
to spatial and temporal fluctuations even in the purely deterministic model 
(Figure [^. 

Figure shows how the number of predators and prey at the two bound¬ 
aries (a: = 0, a: = 20) vary time for the mean held model. We note that 
during the hrst few oscillations both predator and prey populations are close 
to zero. In the corresponding stochastic model, the population can only be 
integer-valued, so a value below 1 in the deterministic model indicates that 
extinction of the population is likely. At later times, we observe regular os- 
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cillations with minima significantly above 0, and conclude that if extinction 
were to occur in the stochastic model, it would most likely happen at early 
times. Figure [^conhr ms these expectations. Figure 8(a) shows that, for the 
stochastic model, the spatial average of the number of predators (or prey) 
may exhibit oscillations similar to those of the deterministic model. How¬ 
ever, extinction of either population can also occur. Figure 8(c) shows that 
if the prey die out hrst, then, necessarily, the predators will also die out as 
well. On the other hand, if the predators die out hrst, then prey numbers 
will blow-up (see Figure [8(e)| . Figures 8(b)[(d)[(f) conhrm that the hybrid 
model can reproduce each of these three qualitatively different scenarios. 
Since statistics for the mean and standard deviation are not meaningful if 
the prey population blows up, we compare the frequency of extinction events 
in the stochastic and hybrid models. 



Stochastic 

Hybrid 100 

Hybrid 50 

Hybrid 25 

Hybrid 10 

Prey 

0.46 

0.45 

0.46 

0.45 

0.43 

Predator 

0.59 

0.57 

0.59 

0.61 

0.56 


Table 2: Shown is the extinction probability for the prey and predator pop¬ 
ulation in the Lotka-Volterra system for the scenario as described in this 
section. This probability is calculated by repeating 256 simulations for each 
of the stochastic model and the hybrid model with thresholds of 10, 25, 50 
and 100, and recording at time f = 80 if a population is extinct or not. 


The results presented in Table suggest that the hybrid model with a 
threshold of 0 = 10 or larger has a similar probability of extinction as the 
stochastic model. We conclude that in this case the hybrid model provides 
a good approximation to the stochastic model. We hnally investigate the 
performance gain obtained by using the hybrid model. 

Table compares the computational time needed to perform 256 simula¬ 
tions of the stochastic and hybrid models proceeding until extinction or time 
t = 100 in all cases. The simulations were performed on a Xeon-2680 16-core 
server with 2.7 GHz and 128GB RAM. Simulations of the stochastic model 
were stopped once predator extinction occurred, as in that case the popula¬ 
tion of prey necessarily blows up. Hence, the advantage in performance of 
the hybrid model is even larger than the numbers indicate. 
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Stochastic 

Hybrid 100 

Hybrid 10 

Simulation time in hours 

40.0 

10.3 

4.8 


Table 3: Performance evaluation of hybrid model. We performed 256 simula¬ 
tions for each of the pure stochastic model as well as the hybrid model with 
a threshold of 10 and 100. As expected, we see a significant performance 
gain when using the hybrid model compared to the stochastic model, and 
the speedup is larger the lower the threshold is, as then we switch earlier to 
the PDE. 
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Figure 6: Simulations of the spatial Lotka-Volterra Model with parameters 
Dn = Dm = 1, a = 2, 6 = 0.1, c = 3, kmax = 101, h = 0.2 and initial values 
N{k, t = 0) = 50, M{k, t = 0) = 5 for fc = 1,..., k^ax- We compare (b)| 
the stochastic model, to[^ (d), the hybrid model with thresholds 0 = 10, 
and (e) (f) with 0 = 25. The hgures on the left show the spatial prohle at 


time t = 4.1 of the number of predators and prey of a single realisation as 
well as the mean of 256 different realisations, whereas the hgures on the right 
show the time evolution of the spatial average of numbers of prey in a single 
realisation and the mean of realisations. The corresponding PDF solution is 
shown in Figure 



















































































Figure 7: Number of prey, left, and predators, right, obtained from the 
deterministic model, equation (24), at the two boundaries, x = 0 and x = 
20, as a function of time. The initial steep rise of the number of prey at 
X = 20 is due to the absence of predators in this region and cannot be 
fully seen here. It is followed by a sharp peak of predators. We see that 
during the first few oscillations the number of individuals can often be close 
to zero, indicating that extinction would be possible at those times in the 
corresponding stochastic model. At later times f > 20 the number of prey 
and predators is fluctuating with minima sufficiently away from zero, making 
extinction at those times less likely. 
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(c) Stochastic Model, Extinction 


6 8 
Time t 

(d) Hybrid Model, Extinction 




Timet 

(e) Stochastic Model, Blow-Up 


Timet 

(f) Hybrid Model, Blow-Up 


Figure 8: Simulations of the spatial Lotka-Volterra Model with parameters 
Dn = Dm = 1, a = 1,6 = 0.1, c = 1, L = 20, kmax = 101. We plot 
the spatial average (iV)^ = Ylk=T -^(^) of fh® number of predators, 
and likewise prey, in a box over time for three realisations of the stochastic 
model, (a)[(c)[(e) and three realisations of the hybrid model, (b)[(d)[(f) with 
a threshold of 0 = 10. These three realisations of each model show the three 
qualitatively different outcomes, namely, oscillatory solutions, extinction of 
both species or extinction of predato^and subsequent blow-up of prey. 




































6. Conclusions 


In this paper we have presented a hybrid algorithm which conples a 
stochastic reaction-diffusion system on a lattice to its associated mean field 
limit, which can be seen as the finite-difference discretisation of a reaction- 
diffusion PDE. Our algorithm preserves mass at the interface between the 
stochastic and deterministic domains, and these domains need neither to 
be static nor connected. Furthermore, for multi-species systems, the cor¬ 
responding stochastic and deterministic domains may differ for individual 
species. We also introduced a normalisation procedure to ensure that the 
stochastic domain contains only integer numbers of particles when the mean 
field equations are evolved or the interface is moved. With the example of 
the spatial Lotka-Volterra system, this paper provides a detailed study of a 
multi-species hybrid system that can accommodate multiple moving inter¬ 
faces. We found that our hybrid algorithm can produce stochastic effects 
that are not present in the corresponding mean field model with the same 
frequency as the stochastic model, while the time taken to perform hybrid 
simulations is much shorter than that for fully stochastic simulations. 

At present, we solve the deterministic equations on the same grid as 
the stochastic model. For computational reasons, one might consider other 
numerical schemes to solve the mean field equations [m ES]. This could 
be particularly important for simulating systems in higher spatial dimen¬ 
sions and geometries with curved boundaries, in such cases, finite difference 
discretisations might not be sufficiently accurate. Another modification to 
determinate the position of our algorithm could be to refine the condition 
used to the interface. For the models studied in this paper, we found that 
counting the number of particles in a box provides a good threshold condition 
for when to use the stochastic, and when to use the mean field, model. This 
is justified as typically stochastic fluctuations scale with the square root of 
the number of particles. However, for some systems one might need to choose 
the domains directly according to the size of the fluctuations. 

Finally, it would be interesting to test our algorithm on reaction-diffusion 
systems involving larger numbers of individuals, such as models of cancer 
growth [miai], angiogenesis [32] or cell polarity [29|. Such considerations 
are beyond the scope of this paper and are therefore postponed for future 
research. 
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Appendix A. Alterations of the Hybrid Algorithm 

Appendix A.l. Large Stochastic Time Step 

A significat problem with continuous diffusion equations is that, in finite 
times, mass can spread arbitrarily far. As a result when the PDE is evolved, 
its solution can leak arbitrarily far into the stochastic domain unless this 
is prevented by the imposittion of an artihcial condition at the interface. 
In our hybrid algorithm, this problem is avoided by using a hnite-difference 
approximation of the PDE with the same lattice size as the stochastic model 
so that on each time step the PDE solution spreads by only one spatial 
compartment. Thus, the PDE solution can, in one time step, enter the 
interface region, but not the stochastic domain. We remark that the infinitely 
fast spread associated with the continuous diffusion equation is not physically 
realistic, it consequently the hnite difference approximation to the diffusion 
PDE is not necessarily less realistic than the PDE itself. 

A problem can arise when the time r to the next event is larger than the 
maximum time step tpde such that the hnite difference scheme converges. 
In this case, we need several iterations of ([^ to evolve the PDE until time r, 
and during this time the PDE solution could leak into the stochastic domain. 
This problem can be avoided by increasing the size of the interface region, or 
one could choose a smaller lattice constant for the PDE regime (compared 
to the stochastic domain). In the examples shown in the present paper, the 
problem of r > tpde rare. When we encountered r > tpde, we chose i 
such that the hnite diherence scheme (|^ converges with time step T We then 
iterated (|^ i times with time step keeping the interface, still consisting of 
one compartment only, hxed. Hence, the interface is treated as a Neumann- 
no-hux boundary during the time step r. We observed that the total changes 
of mass in the interface box is small relative to the total amount of mass 
present, so the error we introduce due to the artihcial Neumann condition is 
negligible. 

A related problem is that when mass moves from the PDE into the in¬ 
terface domain, this, in principle, changes the transition rates of stochastic 
reactions of the interface, and drift terms to the master must be added equa¬ 
tion. However, as discussed above, for the simulations performed in this 
paper the change of mass in the interface compartment is small within one 
time step. Hence, the change in the transition rates one will also be is small 
and can be neglected. 
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Appendix A.2. Small Stochastic Time Step 

If the random Gillespie time step is much smaller than the PDE conver¬ 
gence time step, r ^ tpde, a simple modihcation of our algorithm consists of 
not iterating the mean field domain during every time step, but only after I 
time steps such that for the i Gillespie time steps r*, we have J2i=i T < 'Tpde, 
but the next time step is likely to bring the cumulative time step above tpde- 


Appendix A.3. Deterministic Interface Reactions 

The interface is chosen such that the mean-field equations are sufficiently 
accurate to represent the system at the interface. In section 3.1.1 we have 


described the reactions in the interface compartment in a stochastic way, but 
we could also describe the reactions deterministically by replacing equation 
(13) with 


(Al) 

Then, the transition rates describing stochastic reactions in the interface 
compartment should be set to zero. 


Appendix A.f. Multiple Interfaces for Different Species 

In many situations, the regions where the concentration of one species 
is high, and hence the deterministic PDE description is valid, are different 
for different species. Hence, there is often a requirement to have separate 
interfaces for the different species. We can use the interface condition as in 
section 3.1.1 for each species separately. However, we will now encounter 
regions of space [Lj^,Lif where some species are modelled stochastically 
and others deterministically. Let Ni{k) be stochastic for k = kj^,... ,ki^, 
and n 2 {k) be deterministic in the same interval, where, as before, we identify 
kj^h = Lp, kph = Lp. For simplicity, we focus on a single reaction involving 
only those two species, which in the full stochastic model would be written 
as 

'TNi{k)+pi,N2{k)+p2\Ni{k),N2{k) = R{Nl{k), N2{k)). (A.2) 


As species 2 is modelled deterministically, we replace N 2 {k) by hn 2 {k) and 
obtain 

'TN-j^[k)+pi,hn2{k)+p2\Ni{k),hn2{k) R{Nl{kf hn2{k)f (A.3) 
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Hence, the reaction is still stochastic. However, hn 2 {k) is real valued. The 
situation is thus similar to how reactions were dealt with at a single-species 
interface. The real-valuedness does not cause any problems as we assume 
hn 2 {k) is large in the interface region. If the reaction vector p 2 is negative, 
we have to assume hn 2 {k) S> P 2 , so that a single stochastic reaction cannot 
lead to negative values of n 2 {k). This is easily ensured by choosing the hybrid 
model threshold accordingly high. 

We can regard this reaction as a stochastic source for an otherwise deter¬ 
ministic n 2 {k). The deterministic part of the evolution of n 2 {k) is described 
by 

n2{k,t + At) = n2{k,t) 

+ At {n2{k + l,t) + n2{k - 1, f) - 2n2(fc, t)) + ra {n2{k)'^ , 

kj^ < k < kj^- (A.4) 


The corresponding equations at k = kj^ oi k = kj^ require modihcations of 
the flux term in the same way as the equations in the single species case. 


equation (13). We have denoted reactions which only depend on n 2 {k) as 
r2 {n2{k)). 

The hnite difference equation ([^ was chosen such that it produces exactly 
the mean behaviour of the underlying stochastic model. We now conhrm that 


the mixed model (A.3) and (A.4) still produces the same mean behaviour in 


the limit of large particle numbers. This means we assume Ni{k) oc fl, with 
1, whereas n 2 {k) oc already to justify the use of the deterministic 
equations. We have, by the derivation of (|^, that 


TNi{k)+pl,N2{k)+p^\Ni{k),N2{k) — R{Ni{k), N2{k)) 


= Q (^r{ni{k),n2{k))+ 0 . (A.5) 


But this means we will also have 


R 

as required. 


n2(fc)^ =n(^r {ni{k),n2{k)) + O ^ , (A.6) 
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